###############################################
# NATIVE ADMINISTRATIONS' BUDGETS #############
#
# Replication Material for:
# Continuity or Change? (In)direct Rule in British and French Colonial Africa
# 
# Carl Mueller-Crepon, 2020
# International Organization
#
# File Description:
# Main analysis of native adminsitration bugets,
# local level number of European administrators,
# and the class of chiefs in Nigeria.
#
#
# Called from indrule_analysis_all.R 
#
################################
rm(list = ls())


# INITIALIZE

# Libraries
library(lfe)
library(stargazer)
library(plyr)
library(ggplot2)
library(viridisLite)
library(maptools)
library(rgeos)
library(raster)
library(MASS)
library(lme4)
library(GGally)
library(car)
library(reshape)

# Paths
fig.path = "figures/governance/"
tab.path = "tables/governance/"

# Output type
output.type <- "latex"


##### LaTeX Functionalities
source("scripts/functions/analysis_functions.R")


# Balance table
source("scripts/functions/balance_table.R")

# Stargazer
output.type = "latex"

# Load data
load("data/data_budgetetal.RData")



################################
# BRITISH COLONIES
################################


# DATA #########################

# Interaction terms
data.cross$v33.num_max.cash.crop <- data.cross$v33.num * data.cross$max.cash.crop
data.brit$v33.num_max.cash.crop <- data.brit$v33.num * data.brit$max.cash.crop
data.frnc$v33.num_max.cash.crop <- data.frnc$v33.num * data.frnc$max.cash.crop

# Population density 1880
data.cross$popdens1880.ln <- log((data.cross$pop.1880+1) / data.cross$area)
data.brit$popdens1880.ln <- log((data.brit$pop.1880+1) / data.brit$area)
data.frnc$popdens1880.ln <- log((data.frnc$pop.1880+1) / data.frnc$area)

# Ethnic population density 1880
data.cross$group.popdens.ln <- log((data.cross$group.pop.1880+1) / data.cross$group.area ) 
data.brit$group.popdens.ln <- log((data.brit$group.pop.1880+1) / data.brit$group.area ) 
data.frnc$group.popdens.ln <- log((data.frnc$group.pop.1880+1) / data.frnc$group.area ) 

# Logged distances
log_dist <- function(df){
  for(v in c("river.dist", "coast.dist")){
    df[,paste0(v, ".ln")] <- log(1 + df[,v])
  }
  df
}
data.cross <- log_dist(data.cross)
data.brit <- log_dist(data.brit)
data.frnc <- log_dist(data.frnc)

# For European admin data
data.admin$coast.dist.ln <- log(data.admin$coast.dist + 0.1)
data.admin$river.dist.ln <- log(data.admin$river.dist)
data.admin$popdens1880.ln <- log(data.admin$pop.1880 / data.admin$area)
data.admin$group.popdens.ln <- log((data.admin$group.pop.1880+1) / data.admin$group.area )



# SETUP ##########################

# ... explanatory variables
main.expl <- c("v33.num")

# ....... raw
base.controls <- c("popdens1880.ln", 
                   "group.popdens.ln","coast.dist.ln", "river.dist.ln", "log(1+pop)","log(area)")
base.controls.abs <- c("popdens1880.ln", 
                   "group.popdens.ln","coast.dist.ln", "river.dist.ln")
ethnic.controls <- paste0("v",c(4,5,28),".num")
nature.controls <- c("medianaltitude","medianslope","temperaturemean","evapotranspiration",  
                     "precipitation", "ppetratio","suit", "max.cash.crop" )


# ... fixed effects
fe.panel = " colony.year "
fe.cross = " colony "

# ... SE cluster on spatial unit
se.cluster <- " province "

# ... output specification for stargazer

# ... variables to keep and labels
keep.controls <- c("v33.num", "popdens1880.ln","group.popdens.ln","coast.dist.ln", "river.dist.ln","pop", "area")
keep.labels <- c("Precol. centralization","Pop. density 1880 (log) ","Ethnic pop. density 1880 (log) ", 
                 "Distance to coast (log)", "Distance to river (log)", 
                 "Population (log)", "Area (log)")

keep.controls.abs <- c("v33.num", "popdens1880.ln","group.popdens.ln", "coast.dist.ln", "river.dist.ln")
keep.labels.abs <- c("Precol. centralization","Pop. density 1880 (log) ","Ethnic pop. density 1880 (log) ", 
                     "Distance to coast (log)", "Distance to river (log)")

# ... notes below tables
notes.p <- "$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01."
latex.notes.add <- "Standard errors are clustered on the province-level."
latex.notes.add <- function(width = .8, sample = "The sample includes the colonies of the Gold Coast (Ghana), Nigeria, Nyasaland (Malawi), and Uganda." , cluster = "province-level"){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} OLS models. Standard errors are clustered on the ",cluster,". ", 
         sample,
          " Nature controls consist of median altitude and slope, mean annual temperature, precipitation
         and evapotranspiration, the ratio of the two, agricultural suitability, and soils' suitability for cash crop production. 
         Ethnic controls include the reliance on agriculture and pastoralism, as well as the intensity of agricultural activities. 
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.}")}
latex.notes.add.appx <- function(width = .8, sample = "The sample includes the colonies of the Gold Coast (Ghana), Nigeria, Nyasaland (Malawi), and Uganda.", cluster = "province-level"){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} OLS models. Standard errors are clustered on the ",cluster,". ", 
         sample,
         " Baseline controls consist of the logged 1880 population density of the district and its ethnic groups, 
the logged distance to coast and closest navigable river, and, for per-capita outcomes, the logged district area and population.
          Nature controls are the median altitude and slope, mean annual temperature, precipitation
         and evapotranspiration, the ratio of the two, agricultural suitability, and soils' suitability for cash crop production. 
         Ethnic controls include the reliance on agriculture and pastoralism, as well as the intensity of agricultural activities. 
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.}")}

# Types of revenues and expenditures
rev.types <- c("tax", "fees", "transfer", "other")
exp.types <- c("admin","order","educsoc","agric","works","misc")


#######################################
# DECRIPTIVE PLOTS AND TABLES
# Figures 7 and A7, A8
# Tables 3, 4, A8, A9
#######################################
source("scripts/governance/analysis_plots.R")


#######################################
# FRENCH WEST AFRICA
# Tables A26, A27
#######################################
source("scripts/governance/analysis_wa.R")




#######################################
# EUROPEAN ADMINISTRATORS
# Appendix Table A21
#######################################
stub <- "staff.brit"


# Specification
expl.vars <- c("v33.num")
staff.var <- "I((1e6 * staff/pop+1))"
staff.cluster <- " geo.match"
add.vars <- list(c(base.controls.abs),
                 c(base.controls.abs, nature.controls),
                 c(base.controls.abs, ethnic.controls),
                 c(base.controls.abs, ethnic.controls,  nature.controls))

# Estimation
this.m <- lapply(add.vars, function(av){
  return(felm(make_form(dv = staff.var, expl = c(expl.vars, av),
                        iv = "0", fe = fe.cross, se =  staff.cluster), 
              data = data.admin))
})


# Output

# ... info
add.lines <- list(latex.any.row("Colony FE: ",rep("yes",length(this.m))),
                  latex.any.row("Baseline controls: ", rep(c("yes","yes","yes","yes"),1)),
                  latex.any.row("Nature controls: ", rep(c("no","yes","no","yes"),1)),
                  latex.any.row("Ethnic controls: ", rep(c("no","no","yes","yes"),1)),
                  latex.mean.dv(this.m))

# ... stargazer
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(stargazer(this.m,
                     title="Local-level European Administrators: Nigeria and Uganda",
                     dep.var.caption = "",dep.var.labels.include = FALSE,
                     column.labels = c("European administrators per million"),column.separate = c(4),
                     keep= c("v33.num"),
                     covariate.labels = c("Precol. centralization"),
                     notes.align = "l",label=paste0("tab.",stub),align =T,
                     add.lines = add.lines, 
                     digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),column.sep.width = "4pt",
                     notes = latex.notes.add.appx(sample = "The sample consists of the colonies of Nigeria and Uganda.", cluster = "district-level"), notes.label = "", notes.append = F,
                     type  = output.type),
           fileConn)
close(fileConn)



#######################################
# CHIEF CLASSES
# Appendix Table A28
#######################################
stub <- "chief.class"

# Plot
png(paste0(fig.path, "chief_scatter.png"), 
    height = 3, width = 4, res = 300, unit = "in")
g <- ggplot(data.cross[data.cross$max.chief.class %in% c(0:3),], aes(x = v33.num , y = max.chief.class)) + 
  geom_point() + geom_smooth(method = "lm") +
  ylab("Class of chiefs") + xlab("Precolonial centralization")
print(g)
dev.off()


# Specification
expl.vars <- c("v33.num")
add.vars <- list(c(base.controls.abs),
                 c(base.controls.abs, nature.controls),
                 c(base.controls.abs, ethnic.controls,  nature.controls))

# Estimation
this.m <- lapply(add.vars, function(av){
  return(felm(make_form(dv = "max.chief.class", expl = c(expl.vars, av),
                        iv = "0", fe = fe.cross, se =  " geo.match"), 
              data = data.cross))
})


# Output

# ... info
add.lines <- list(latex.any.row("Colony FE: ",rep("yes",length(this.m))),
                  latex.any.row("Baseline controls: ", rep(c("yes","yes","yes"),1)),
                  latex.any.row("Nature controls: ", rep(c("no","yes","yes"),1)),
                  latex.any.row("Ethnic controls: ", rep(c("no","no","yes"),1)),
                  latex.mean.dv(this.m))

# ... stargazer
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(stargazer(this.m,
                     title="Highest class of chief in district: Nigeria (1924/1929)",
                     dep.var.caption = "",dep.var.labels.include = FALSE,
                     column.labels = c("Highest class of chief (1-3)"),column.separate = c(6),
                     keep= c("v33.num"), 
                     covariate.labels = c("Precol. centralization"),
                     notes.align = "l",label=paste0("tab.",stub),align =T,
                     add.lines = add.lines, 
                     digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),column.sep.width = "5pt",
                     notes = latex.notes.add.appx(sample = "The sample consists of the colony of Nigeria.", cluster = "district-level"), notes.label = "", notes.append = F,
                     type  = output.type),
           fileConn)
close(fileConn)


###################################
# REVENUES AND EXPENDITURES: Main Table
# Table 6, main paper
###################################

stub <- "revexp.comb"

# Estimate
this.m <- list(felm(make_form(dv = "log(total.rev +1)", expl = c(main.expl, base.controls.abs,  ethnic.controls,  nature.controls),
                    iv = "0", fe = fe.cross, se =  se.cluster), 
          data = data.cross),
     felm(make_form(dv = "log((total.rev + 1)/pop)", expl = c(main.expl, base.controls,  ethnic.controls,  nature.controls),
                    iv = "0", fe = fe.cross, se =  se.cluster), 
          data = data.cross),
     felm(make_form(dv = "log(total.exp +1)", expl = c(main.expl, base.controls.abs,  ethnic.controls,  nature.controls),
                    iv = "0", fe = fe.cross, se =  se.cluster), 
          data = data.cross),
     felm(make_form(dv = "log((total.exp + 1)/pop)", expl = c(main.expl, base.controls,  ethnic.controls,  nature.controls),
                    iv = "0", fe = fe.cross, se =  se.cluster), 
          data = data.cross))



# Output

# ... info
add.lines <- list(latex.any.row("Colony FE: ",rep("yes",length(this.m))),
                  latex.any.row("Nature controls: ",rep("yes",length(this.m))),
                  latex.any.row("Ethnic controls: ",rep("yes",length(this.m))),
                  latex.mean.dv(this.m))
add.sec.row.head <- "} \\\\ & \\multicolumn{1}{c}{total} & \\multicolumn{1}{c}{per capita}  & \\multicolumn{1}{c}{total} & \\multicolumn{1}{c}{per capita"

# Save
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(stargazer(this.m,
                     title="Native treasuries under British rule: per-capita revenues and expenditures (logged 2016 $\\pounds$)",
                     dep.var.caption = "",dep.var.labels.include = FALSE,
                     column.labels = collab_w_ruler(column.labels = c("Revenues", "Expenditures"), 
                                                    column.separate = c(2,2), trim = 0, add.below = rep(c("total","per capita"), 2)),
                     column.separate = c(2,2),
                     keep= keep.controls, covariate.labels = keep.labels,
                     notes.align = "l",label=paste0("tab.",stub),align =T,
                     add.lines = add.lines, 
                     digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),column.sep.width = "-0pt",
                     notes = latex.notes.add(.95), notes.label = "", notes.append = F,
                     type  = output.type),
           fileConn)
close(fileConn)


###################################
# REVENUES BY TYPE
# Appendix Table A24
###################################

stub <- "rev.type"
rev.types <- c("tax", "fees", "transfer", "other")
rev.titles <- c("Taxes", "Fees \\& Fines", "Transfers", "Other")

# Specification
dep.var <- paste0("I(log((",rev.types,"+1)/pop))")

spec <- unlist(lapply(dep.var, function(dv){
  make_form(dv = dv, expl = c(main.expl, base.controls,  ethnic.controls,  nature.controls),
            iv = "0", fe = fe.cross, se =  se.cluster)}))

# Estimation
this.m <- lapply(spec, function(f){
  return(felm(f, data = data.cross))
})
# Save for plotting
m.revtypes <- this.m

# Output

# ... info
add.lines <- list(latex.any.row("Colony FE: ",rep("yes",length(this.m))),
                  latex.any.row("Baseline controls: ",rep("yes",length(this.m))),
                  latex.any.row("Nature controls: ",rep("yes",length(this.m))),
                  latex.any.row("Ethnic controls: ",rep("yes",length(this.m))),
                  latex.mean.dv(this.m))

# ... stargazer
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(stargazer(this.m,
                     title="Native treasury revenues per capita by type (2016 $\\pounds$)",
                     dep.var.caption = "Revenues/capita (log)",dep.var.labels.include = FALSE,
                     column.labels = c("Taxes","Fees \\& fines", "Transfers","Other"),
                     keep= keep.controls[1], covariate.labels = keep.labels[1],
                     notes.align = "l",label=paste0("tab.",stub),align =T,
                     add.lines = add.lines, 
                     digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.add.appx(.86), notes.label = "", notes.append = F,
                     type  = output.type),
           fileConn)
close(fileConn)


###################################
# EXPENDITURES BY TYPE
# Appendix Table 25
###################################

stub <- "exp.type"
exp.types <- c("admin","order","educsoc","agric","works","misc")
exp.titles <- c("Admin.","Order","Educ. \\& Health","Agric.","Works","Other")

# Specification
dep.var <- paste0("I(log((",exp.types,"+1)/pop))")

spec <- unlist(lapply(dep.var, function(dv){
  make_form(dv = dv, expl = c(main.expl, base.controls,  ethnic.controls,  nature.controls),
            iv = "0", fe = fe.cross, se =  se.cluster)
}))

# Estimation
this.m <- lapply(spec, function(form.str){
  return(felm(as.formula(form.str), data = data.cross))
})
# Save for plotting
m.extypes <- this.m

# Output

# ... info
add.lines <- list(latex.any.row("Colony FE: ",rep("yes",length(this.m))),
                  latex.any.row("Baseline controls: ",rep("yes",length(this.m))),
                  latex.any.row("Nature controls: ",rep("yes",length(this.m))),
                  latex.any.row("Ethnic controls: ",rep("yes",length(this.m))),
                  latex.mean.dv(this.m))

# ... stargazer
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(stargazer(this.m,
                     title="Native treasury expenditures per capita by type (2016 $\\pounds$)",
                     dep.var.caption = "Expenditures/capita (log)",dep.var.labels.include = FALSE,
                     column.labels = exp.titles,
                     keep= keep.controls[1], covariate.labels = keep.labels[1],
                     notes.align = "l",label=paste0("tab.",stub),align =T,
                     add.lines = add.lines, 
                     digits = 2, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),column.sep.width = "-7pt",
                     notes = latex.notes.add.appx(.95), notes.label = "", notes.append = F,
                     type  = output.type),
           fileConn)
close(fileConn)



#######################################
# ROBUSTNESS CHECKS
#######################################

#######################################
# Main table with:
# 1. Desease Environment
# 2. Weight by colony 
# 3. Drop outliers 
# 4. Panel regression
# 5. Hierarchical model
#
# Appendix Table A22
#######################################

# ... Stub
stub <- "revexp.robust1"

# ... dependent variable
dep.var <- c("log((total.rev + 1)/pop)")

# ... Make empty list
this.m <- list()

# ... Additional controls: desease environment
desease.vars <- c("malaria.tsuit.max","tsetse.max")
this.m <- c(this.m, 
            list(felm(make_form(dv = dep.var, expl = c(main.expl,base.controls,   desease.vars,  ethnic.controls,  nature.controls),
                                iv = "0", fe = fe.cross, se =  se.cluster), 
                      data = data.cross)))


# ... Weight by colony
form <- make_form(dv = dep.var, expl = c(main.expl,base.controls,  ethnic.controls,  nature.controls),
                  iv = "0", fe = fe.cross, se =  se.cluster)
this.data <- na.omit(data.cross[, all.vars(form)] )
this.data$weights <- gen_weights(this.data$colony)
this.m <- c(this.m, 
            list(felm(form, data = this.data, weights = this.data$weights)))


# ... Drop outliers
this.m <- c(this.m, 
            list(felm(make_form(dv = dep.var, expl = c(main.expl,base.controls,  ethnic.controls,  nature.controls),
                                iv = "0", fe = fe.cross, se =  se.cluster), 
                      data = data.cross[data.cross[,"total.rev"] > quantile(data.cross[,"total.rev"],.025, na.rm = T) &
                                          data.cross[,"total.rev"] < quantile(data.cross[,"total.rev"], .975, na.rm = T),])))

# ... Panel Regression, district weights
form <- make_form(dv = dep.var, expl = c(main.expl,base.controls,  ethnic.controls,  nature.controls),
                  iv = "0", fe = "factor(colony):factor(year)", se =  se.cluster)
this.data <- na.omit(data.brit[, c("geo.match",all.vars(form))] )
this.data$weights <- gen_weights(this.data$geo.match)
this.m <- c(this.m, 
            list(felm(form, data = this.data, weights = this.data$weights)))

# ... Hierarchical model

# ... Specification
random.effect <- " + (1| geo.match) + (1|colony)"
fixed.effect <- "+ factor(colony.year)"

form.str <- paste0(dep.var, "~", paste(c(main.expl,
                   base.controls, ethnic.controls,nature.controls), collapse = " + "),
                   fixed.effect, random.effect)
this.m <- c(this.m, list(lmer(as.formula(form.str), data = data.brit)))



# Output

# ... info
add.lines <- list(latex.any.row("Fixed effect: ",c("colony","colony","colony","col.-year","col.-year")),
                  # latex.any.row("Random effects: ",c("--","--","--","district \\&","--","--")),
                  # latex.any.row("",c("--","--","--","colony","--","--")),
                  # latex.any.row("Weigths: ",c("--","colony","district","--","--","")),
                  # latex.any.row("Disease controls: ",c("yes","no","no","no","no","no")),
                  latex.any.row("Baseline controls: ",rep("yes",length(this.m))),
                  latex.any.row("Nature controls: ",rep("yes",length(this.m))),
                  latex.any.row("Ethnic controls: ",rep("yes",length(this.m))),
                  latex.mean.dv(this.m))

# ... stargazer
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(stargazer(this.m,
                     title="Per-capita revenues: Robustness checks",
                     dep.var.caption = "Revenues p.c. (log)",dep.var.labels.include = FALSE,
                     column.labels = c("Desease","Col.-weight","No outlier", "Wght. panel","HLM"),
                     keep= keep.controls[1] , covariate.labels = keep.labels[1],
                     notes.align = "l",label=paste0("tab.",stub),align =T,
                     add.lines = add.lines, model.names = F,
                     digits = 2, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"), column.sep.width = "-0pt",
                     notes = "\\parbox[t]{0.98\\textwidth}{\\textit{Notes:} OLS models in 1--4, 
                      hierarchical linear model in 5. The sample includes the colonies of the Gold Coast (Ghana), Nigeria, Nyasaland (Malawi), and Uganda. 
                     Standard errors are clustered on the province-level.
Baseline controls are the logged 1880 population density of the district and its ethnic groups, 
the logged distance to coast and closest navigable river, and the logged district area and population.
                     Nature controls include the median altitude and slope, mean annual temperature, precipitation
                     and evapotranspiration, the ratio of the two, agricultural suitability, and soils' 
                     suitability for cash crop production. 
                     Ethnic controls are the reliance on agriculture and pastoralism, 
                     as well as the intensity of agricultural activities. 
                     Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.}", notes.label = "", notes.append = F,
                     type  = output.type),
           fileConn)
close(fileConn)

#######################################
# Alternative indirect rule specifications
# 1. M & P Coding
# 2. Ruler 1885
# 3. Max class of chief
#
# Appendix Table A23
########################################

# ... Stub
stub <- "revexp.altind"

# ... Specifications
spec <- unlist(lapply(c("v33.num.mp","has.ruler.1885","max.chief.class"), function(iv){
  c(make_form(dv = "log((total.rev + 1))", expl = c(iv,base.controls.abs,  ethnic.controls,  nature.controls),
              iv = "0", fe = fe.cross, se =  se.cluster),
    make_form(dv = "log((total.rev + 1)/pop)", expl = c(iv,base.controls,  ethnic.controls,  nature.controls),
              iv = "0", fe = fe.cross, se =  se.cluster))
}))


# ... Estimation
this.m <- lapply(spec, function(form){
  return(felm(form, data = data.cross))
})



# Output 

# ... info
add.lines <- list(latex.any.row("Colony FE: ",rep("yes",length(this.m))),
                  latex.any.row("Baseline controls: ",rep(c("yes","yes"),length(this.m)/2)),
                  latex.any.row("Nature controls: ",rep(c("yes","yes"),length(this.m)/2)),
                  latex.any.row("Ethnic controls: ",rep(c("yes","yes"),length(this.m)/2)),
                  latex.mean.dv(this.m))

# ... stargazer
fileConn<-file(paste0(tab.path, stub, ".tex"))
writeLines(stargazer(this.m,
                     title="Revenues  (2016 $\\pounds$): Alternative specifications",
                     dep.var.caption = "Revenues (log):",dep.var.labels.include = FALSE,
                     column.labels = c("Total","Per capita","Total","Per capita","Total", "Per capita"),
                     keep= c("v33.num.mp", "has.ruler.1885","max.chief.class"),
                     order= c("v33.num.mp", "has.ruler.1885","max.chief.class"),
                     covariate.labels = c("Precol. centr. (M\\&P)", "Capital 1885", "Chief class"),
                     notes.align = "l",label=paste0("tab.",stub),align =T,
                     add.lines = add.lines, column.sep.width = "-10pt",
                     digits = 2, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.add.appx(.95), notes.label = "", notes.append = F,
                     type  = output.type),
           fileConn)
close(fileConn)


######################################
# Colony-Jackknife
# Appendix Figure A15
######################################

stub <- "revexp.coljack"

# ... colonies / countries
colonies <- unique(data.cross$colony)
colonies <- c("Baseline", colonies[order(colonies)])

# ... estimation
this.m <- unlist(lapply(colonies, function(c){
  list(felm(make_form(dv = "log(total.rev +1)", expl = c(main.expl, base.controls.abs),
                      iv = "0", fe = " 0 ", se =  "geo.match"), 
            data = data.cross[data.cross$colony != c,]),
       felm(make_form(dv = "log(total.rev + 1)", expl = c(main.expl, base.controls.abs,  nature.controls),
                      iv = "0", fe = " 0 ", se =  "geo.match"), 
            data = data.cross[data.cross$colony != c,]),
       felm(make_form(dv = "log(total.rev +1)", expl = c(main.expl, base.controls.abs,  ethnic.controls,  nature.controls),
                      iv = "0", fe = " 0 ", se =  "geo.match"), 
            data = data.cross[data.cross$colony != c,]),
       felm(make_form(dv = "log(total.rev +1)", expl = c(main.expl, base.controls.abs),
                      iv = "0", fe = fe.cross, se =  "geo.match"), 
            data = data.cross[data.cross$colony != c,]),
       felm(make_form(dv = "log(total.rev + 1)", expl = c(main.expl, base.controls.abs,  nature.controls),
                      iv = "0", fe = fe.cross, se =  "geo.match"), 
            data = data.cross[data.cross$colony != c,]),
       felm(make_form(dv = "log(total.rev +1)", expl = c(main.expl, base.controls.abs,  ethnic.controls,  nature.controls),
                      iv = "0", fe = fe.cross, se =  "geo.match"), 
            data = data.cross[data.cross$colony != c,]))
}), recursive = F)

# ... get coefficients
plot.df <- data.frame(do.call(rbind,lapply(this.m, function(m){
  x <- c(coef = m$coefficients["v33.num",1], se = diag(m$clustervcv)["v33.num"]^.5)
  x <- matrix(x, nrow = 1, byrow = T)
  colnames(x) <- c("coef","se")
  x
})))

# ... Add names
plot.df$name <- rep(colonies, each = 6)
plot.df$spec <- rep(c("(1) Baseline (abs)", "(2) + nature controls (abs)","(3) + nature +\nethnic controls (abs)"), 
                    length(colonies) * 2)
plot.df$FE <- rep(rep(c("no", "yes"), each = 3),  length(colonies))

# ... change names to current country names
plot.df$name[grep("Gold Coast", plot.df$name)] <- "Ghana"
plot.df$name[grep("Nyasaland", plot.df$name)] <- "Malawi"


# ... sort by name
# plot.df <- plot.df[order(plot.df$name, plot.df$spec),]
plot.df$pos <- rep(c(length(colonies):1), each = 6) + rep(rep(c(.1, -.1), each = 3),  length(colonies))



# ... plot
png(paste0(fig.path,stub,".png"), width = 7, height = 3, unit = "in", res = 600)
g <- ggplot(plot.df, aes(y = pos, x = coef, color = FE)) + geom_point() + 
  geom_errorbarh(aes(xmin = coef - 1.96*se, 
                     xmax = coef + 1.96*se), 
                 height=.1) +
  geom_vline(xintercept = 0, color = "darkgrey", lty = 2) +  
  scale_y_continuous(breaks=rev(c(length(colonies):1)), labels=rev(unique(plot.df$name))) + 
  scale_color_discrete(name = "Colony fixed effects:") +
  facet_wrap(~ spec, nrow = 1) + xlab("Coefficient of pre-colonial centralization") + ylab(NULL) +
  theme(legend.position = "top")
print(g)
dev.off()

